The mismatch between experimental and computational fluid dynamics analyses for magnetic surface microrollers

Magnetically actuated Janus surface microrollers are promising microrobotic platform with numerous potential biomedical engineering applications. While the locomotion models based on a "rotating sphere on a nearby wall" can be adapted to surface microrollers, real-world dynamics may differ from the proposed theories/simulations. In this study, we examine the locomotion efficiency of surface microrollers with diameters of 5, 10, 25, and 50 µm and demonstrate that computational fluid dynamics simulations cannot accurately capture locomotion characteristics for different sizes of microrollers. Specifically, we observe a significant mismatch between lift forces predicted by simulations and opposite balancing forces, particularly for smaller microrollers. We propose the existence of an unaccounted force component in the direction of lift, which is not included in the computational fluid dynamics simulations. Overall, our findings provide a deeper understanding of the physical mechanisms underlying surface microroller locomotion and have important implications for future applications in biomedical engineering.

www.nature.com/scientificreports/ diameters, 5, 10, 25, and 50 µm. We analyzed their translational speeds (V) in a static and semi-infinite environment for increasing rotational speeds (Ω). Our investigation revealed a non-linear relationship between microroller size and locomotion efficiency, defined as the microroller's ability to convert rotational motion into translational motion. The microroller with a diameter of 50 µm demonstrated the highest efficiency among the different sizes tested, while the 10 µm microroller was found to be the least efficient. The overall ranking of the microrollers from most to least efficient was as follows: 50 µm > 5 µm > 25 µm > 10 µm. To understand this relationship, we conducted computational fluid dynamics (CFD) analyses. However, the analyses failed to explain the irregular relationship observed in the experimental results. In fact, the CFD analyses significantly overestimated the forces involved, with F G dominating over F L and F rep . This discrepancy indicated that the force balance on the z-axis was not satisfied. Overall, the main aim of this study is to analyze this discrepancy and discuss the reasons behind it.

Results
The basic speed analyses of microrollers with different diameters. We characterized the speeds of magnetic surface microrollers for different sizes (Fig. 2). The fabrication of surface microrollers is summarized in Fig. 2a. The silica particles with different sizes were monolayered on a glass substrate, and then the particles were sputter-coated with Ni and Au. To obtain ≈5 μm microroller, 4 μm silica particles are sputter coated with 1000 nm Ni and 50 nm Au; the ≈10 μm microroller also has 1000 nm Ni and 50 nm Au on top of 10 μm particles (Fig. 2b). The ≈25 and ≈50 μm microrollers contain 1800 nm Ni and 50 nm Au on top of 25 and 50 μm silica particles, respectively (Fig. 2b). Then, we characterized their steady-state translational speeds for increasing rotational frequencies until f = 80 Hz since it was the highest actuation frequency for 50 μm ( Figure S1); therefore, f = 80 Hz was selected for the highest limit for all groups for a concise comparison. The other microrollers did not also step-out until f = 80 Hz as their translational speed never decreased with increasing rotation frequency (Fig. 2c). The 50 μm microrollers performed best with the highest body length per second (Fig. 2c), and 5 μm microrollers had the second ranking. The performances of the 10 and 25 μm microrollers were close to each other; the 25 μm had slightly higher relative velocity (Fig. 2c). The relative velocity of the microrollers, their translational speed over their diameter, V/2a, was determined by the lubrication distance, δ ′ = δ/a . The microrollers would have the same relative velocity at a fixed δ ′ for a specific rotation frequency; however, δ ′ is deterministically reached for a specific frequency and microroller size depending on the forces on z-axis. Thus, the relative speed differences were due to different δ ′ values for different sizes; smaller δ ′ results in more efficient microroller locomotion. Therefore, we could say that the 50 μm microroller had the smallest δ ′ , while 10 μm had the highest (Fig. 2d), and the overall ranking goes as δ ′ 10µm > δ ′ 25µm > δ ′ 5µm > δ ′ 50µm . The value of δ ′ depends on the force balance of the microroller on the z-axis (Fig. 1).
The governing forces of microrollers at different sizes in the CFD environment. We performed CFD simulations to capture the forces with hydrodynamic origins, such as F P , F D, and F L, for different microroller sizes. We performed CFD simulations in 3D with two different configurations (Fig. 3a), pure rotation and translation, where the natural motion of the microroller is a combination of the two. We can investigate the system by decomposing these two cases due to the linearity of motion and stress Eqs. 15 . We quantified the forces acting on the microroller in both x-and z-axes. The forces on the x-axis in the pure rotation and translation case correspond to F P and F D , respectively 2,15 . The forces in the z-axis give the F L , while the translation and rotation have different contributions to the total F L . We performed the simulations for all microroller sizes where the environments were adjusted to 40a × 40a × 40a to avoid any confinement effects 2 . We rotated all microrollers for f = 30, 50, 75 and 100 Hz and V = 30, 50, 75 and 100 body length per second for different δ ′ values ( Fig. 3b-d) to be able to derive the expressions for the governing forces. The results revealed that the relative F P and F D values overlapped for all sizes (Fig. 3b,c), implying that these two have no size scaling effect. Furthermore, we used a finer mesh for the simulations than the previous work 2 that was validated from the work from Goldman et al. 15 to ensure more reliable results. Therefore, we obtained new expressions for F P and F D with modified correction factors (R 2 > 0.99 for both expressions): Figure 1. The basic force balance on a microroller near a wall. A microroller, rotating with an angular velocity of Ω, with a distance of δ from the nearby wall, creates a propulsion force, F P , balanced out the drag force, F D , on the x-axis. The lift force, F L , and repulsion force, F rep , are balanced by the gravitational over-buoyancy force, F G . The origin of the forces was marked with the color code. www.nature.com/scientificreports/ On the other hand, the lift forces were negligible for the size below 25 μm, and very prominent for 50 μm (Fig. 3d). The other important finding is that the translational motion did not cause any F L , demonstrating the rotational force is the mere contributor to the lift force. We also derived an expression for F L with the new correction factors using the result of 50 μm since CFD loses the sensitivity in F L for smaller sizes (R 2 > 0.99): We found out that F P and F D had no scaling effect, while F L started appearing with bigger microroller sizes. When we quantified rotational (Re Ω ) and translational (Re V ) Reynolds numbers based on the experiments in Fig. 2c (Fig. 3e,f), we observed that the Re Ω for 25 and 50 μm reaches around 0.1-0.4, indicating that inertial forces played some role on the locomotion of such microrollers, while the contribution of translational locomotion was still small. This also agreed well with the F L quantification (Fig. 3d), since the nature of F L is inertial 26 . Moreover, it is worth highlighting that the expressions for propulsion force (F P ) and drag force (F D ) put forth for the low Reynolds number regime 15 remained valid even when some degree of inertia was introduced into the system for the range of parameters studied (Fig. 3b,c,e,f). To better understand the experimental behavior, we now turn our attention to the effects of other forces.
Generic force relations between the governing forces. As demonstrated in the previous section, the F L only becomes noticeable in the CFD environment for microrollers with diameters greater than 25 μm, with the most significant effect observed for 50 μm microrollers. This suggests that larger microrollers may generate greater F L and, therefore may have a higher δ ′ . However, other forces in the force balance may counteract the effects of F L and alter this relationship. In this section, we explored the relationship between forces acting in the z-axis for microrollers of various sizes and frequencies. This section investigates the relationship between the forces in the z-axis for different sizes and frequencies. We calculated the F G using the experimentally calculated with different diameters were prepared on glass substrates. Then magnetic material (Ni) and passivating layer (Au) were sputter-coated on the particles. After sputter coating, the magnetic film was programmed in the out-of-plane direction by applying a 1.7 T uniform magnetic field while the Janus particles were still on the glass substrate. The schematic was prepared with Sketchup 2023 software, https:// www. sketc hup. com. (b) Light microscopy images from surface microrollers with different sizes ≈ 5, 10, 25, and 50 μm diameters were used in the study. The scale bar is 20 μm. (c) The average relative translational speeds of the microrollers depend on increasing rotation frequencies. The most efficient microroller depending on the body length per second was 50 μm and then 5 μm. The least efficient one was 10 μm. (d) A summary of the experimental results. The most efficient microroller was the biggest one, but there was a non-monotonic relation between diameters and speeds, which has to do with their relative lubrication distance. www.nature.com/scientificreports/ densities of the microrollers and interpolated the sizes between them. F P and F L were also calculated for all sizes using Eqs. (2) and (4) for f = 30 and 80 Hz for different δ ′ . F rep was approximated to be the difference between electrostatic repulsion forces (F es ) and van der Waals attraction forces (F vdW ) such that F rep = F es + F vdW (Sup. Note 1) for the same δ . The computed forces in the z-axis were normalized to the corresponding F P values for a concise comparison between different sizes. The results are summarized in Fig. 4. In this section, we demonstrated the results for the diameter of the microroller vs the force relation (e.g., F G /F P ) for different δ ′ . Since the experimental value of δ ′ is unknown and difficult to measure 28 , we presented broad range of δ ′ . Also, it is known that δ ′ increases with increasing f 2,28 ; therefore, the higher δ ′ represents high relative f region and vice versa. As expected, the gravitational forces increased with the microroller diameter (Fig. 4a). They were more dominant at f = 30 Hz due to less F P generated and even went up to roughly four times than F P for 50 μm at δ ′ = 0.1. The effect of F G disappeared remarkably for smaller diameters, especially at small lubrication distances. F G seemed less important at higher frequencies, such as f = 80 Hz (Fig. 4a). F L was also not prominent for smaller diameters and frequencies (Fig. 4b). The magnitude of F L is insignificant compared to F P almost for all sizes, and it went a maximum of 20% of the F P at 50 μm, δ ′ = 0.1 (Fig. 4b). One should note that the y-axis of the graphs in Fig. 4a,b are significantly different, implying that F G is a more dominant force than F L . It was also revealed that the components of F rep was only prevalent when the size is less than 3-4 μm and the microroller was very close to the wall (Fig. 4c,d). As F rep is a surface-related chemical force, we expected it to become more significant as the surface-to-volume ratio increased. Due to the size range of our microrollers, we did not observe a significant contribution from F rep and therefore did not include it in our further discussions.
As in our study, the forces on the z-axis must balance each other on steady-state locomotion; however, that was not the case for our system (Fig. 1, Fig. 4e). As a result, the magnitude F L is much smaller than the F G for all conditions (Fig. 4e), reaching a maximum of ≈10% for f = 80 Hz at higher microroller diameters. Since F G is a permanent force and thus cannot change, the F L generated by the microrollers seemed to be missing in CFD simulations. There is also a discrepancy in theoretical analysis, the division of F L over F G yields: www.nature.com/scientificreports/ When the diameter (2a) of the microroller increases at constant , the δ ′ also tries to increase itself to keep the equation in balance. However, this was not the case for experiments, and there was an irregular relation between size and the translational speed of the microrollers. Overall, it is evident that there is a missing force on the z-axis, which CFD analyses could not capture.
The missing/balancing force on the z-axis. We propose a new force balance for our system, demonstrated in Fig. 5a. A new force must fill the missing force, F b , to balance out the z-axis. We compared F b with F P for different rotation frequencies to get an initial glance. The analyses revealed that F b is a significant force, having similar and even greater values than F P (Fig. 5b). Then we compared the F b with F G ; the analysis demonstrated that the missing/balancing force, F b , is so significant that it is almost entirely missing for 5 μm microroller (Fig. 5c). The F b followed a regular pattern depending on the size of the microroller when compared to F G , meaning that it was the biggest for 5 μm and smallest for 50 μm, even though the ratios were significantly high for all groups (Fig. 5c). Compared to the F L , F b is much greater than F L , showing that the magnitudes of F L are very small in CFD (Fig. 5d).

Discussions
In this work, we found out that the forces on the z-axis were not balancing each other out in the theoretical analyses. The missing force was in the direction of the lift force, F L , which had a similar order with F G . The possible reasons why CFD could not fully capture the experimental characteristics are as follows: 1. The accuracy of CFD simulations can be limited by assumptions such as no-slip boundary conditions 29,30 , which may not be optimal for mesoscale cases, as found here 30 . This approximation can lead to under-/overestimation of forces and the omission of important features such as atomic and surface forces, wetting, and physio-chemical parameters 29,30 . Therefore, for example, As a result, the quantified lift force in simulations may not reflect the actual behavior of microrollers. For example, previous experimental work has shown that 10 μm microroller lifts itself with increasing f , in other words, δ ′ has increased with increasing f, at low Reynolds number regime where lift forces do not exist in CFD simulations 28 . 2. Although CFD simulations are useful for predicting the behavior of microrollers, they may oversimplify the geometry of the objects and overlook important experimental properties. For example, the surface roughness of the nearby wall and the Janus structure of microrollers are not always considered in these simulations 15,29 . For example, recent study by Rinehart et al. has demonstrated that a significant lift force arises from a rotating Janus particle in the Stokes regime due to slip inhomogeneity in the particle body 31 . Other than roughness or body inhomogeneity, slight shape anisotropies in the body of the microrollers due to the thin film deposition may be a contributing factor. To study the effect of anisotropy, we briefly modeled anisotropic microrollers and observed that the contribution of experimental anisotropy is almost negligible (Sup. Note 2). Also, wet friction caused by the surrounding fluid to the microrollers could be the reason for such behavior 32 . 3. It was shown in the study that inertial forces arose for 25 and 50 μm, evidenced by Reynolds numbers (Fig. 3e,f), which could also be the source of non-linearities presented here. 4. The study approximated the repulsion forces as a combination of electrostatic repulsion force and van der Waals forces. However, it is important to note that other fundamental forces are present in the system, such as solvation forces, including hydrophobic/hydrophilic interactions or other steric effects, which were not taken into account. The origin of hydration or solvation forces differs from that of electrostatic or van der Waals forces, as they arise from the liquid medium and the physical and chemical properties of the interacting surfaces. The effect of these forces is less understood 27 . To accurately determine these forces, rigorous experimentation is required, and they cannot be modeled without determining experimental constants. Therefore, the presence of these forces could also contribute to the observed discrepancy, particularly at smaller lubrication distances.
Ideally, performing a dynamic-multiphysics simulation, including chemical and physical attraction forces and solid-state mechanics, could improve the results; nevertheless, it should be noted that it is exceptionally challenging to make a physical model capturing all features of the surface microrollers. Overall, the results presented here provided an improved understanding of microroller research from a practical perspective with their potential in drug delivery and lab-on-a-chip applications.

Materials and methods
Fabrication and actuation of surface microrollers. Magnetically actuated, spherical Janus microrollers were fabricated by sequentially sputtering Ni and Au nanofilms on pre-dried monolayer of silica (SiO 2 ) particles of 4, 10, 25 and 50 μm diameter (Microparticles GmbH for 4, 10 and 25 μm, Corpuscular Inc.) using benchtop sputter coating system (Leica EM ACE600, Leica Microsystems). 1000 nm Ni was used for 4 and 10 μm, 1800 nm was Ni used for 25 and 50 μm particles, all microparticles had 50 nm Au as passivating layers. After sputtering, magnetization direction of the Janus microparticles was oriented towards out of the cap by applying a 1.8 T uniform magnetic field in a vibrating-sample magnetometer (VSM; MicroSense, Lowell, MA). Ni and Au nanofilm-sputtered microparticles were then released from the substrate via sonication in ethanol. www.nature.com/scientificreports/ The microrollers were intensely washed and finally dispersed in PBS 1 × . The microrollers were actuated using a custom-made five-coiled electromagnetic coil system placed on a microscope (Zeiss Axio Observer A1, Carl Zeiss). The microrollers were actuated using uniform rotating magnetic fields with an amplitude of 10 mT that enables surface rolling and steering control. All the experiments were performed in PBS 1 × with no confinements from any direction. The steady state translational speeds of the microrollers were analyzed using an inhouse MATLAB tracking code.
Computational fluid dynamics (CFD) analyses. COMSOL Multiphysics 6.0 Simulation Software (COMSOL, Inc.) was used to simulate microrollers and calculate the forces acting on the microrollers in all conditions, using laminar flow interface physics by solving the Navier-Stokes equations. We used previously validated simulations with increased resolution 2 . We defined two different mesh regions, the first one is around the microroller and the second one covered the remaining spaces in the simulation environment. 3a × 3a × 3a, where minimum element size, maximum element size, maximum element growth rate, curvature factor, and the resolution of narrow regions were defined in COMSOL as a/50, a/5, 1.2, 0.05, and 1.5 respectively for all sizes, where a is the radius of the microroller. For the second mesh, "extremely fine" built-in mesh settings were used in 40a × 40a × 40a dimensions. The density and the dynamic viscosity of the fluid were taken as ρ=1000 kg/m 3 and µ = 1 cP. The boundaries other than microrollers' were defined as no-slip boundaries. We either rotated or translated the spherical microrollers with diameters, 2a, of 5, 10, 25 and 50 µm to derivate the expressions for F P , F D and F L . The hydrodynamic forces acting on the spheres were quantified in each case according to: where σ is the stress tensor, n is the outward pointing normal vector, and S is the particle surface. The total acting force resulted in + x direction with clockwise rotation of the body, which is propulsion force, F P and + z direction is the contribution of F L from rotation 2 . The total acting force resulted in -x direction with the translation of the microrollers in + x direction is drag force, F D , and + y direction is the contribution of F L from translation. The lift forces on the anisotropic shapes were calculated by taking the average of the forces over one rotation cycle.  The relation between gravitational force (F G ) and propulsion force (F P ) for 30 and 80 Hz for different δ ′ . Gravitational forces start to become more prominent at greater microroller diameters. (b) The relation between the lift force (F L ) and the propulsion force (F P ). At smaller diameters, lift forces were less pronounced. At higher diameters, at δ ′ =0.1, only 20% of the propulsion force was generated, unlike the gravitational forces. (c) The relation between electrostatic repulsion forces and propulsion forces. Electrostatic forces were only prominent when the microroller was very close to the wall, δ ′ =0.005, and became invisible at higher microroller diameters and lubrication distances; therefore, it is assumed to be negligible. (d) The relation between van der Waals attraction forces and propulsion forces. The magnitude of the forces were negligible in the size regime of 5-50 μm; therefore, their effect is assumed to be negligible. (e) The relation between the lift force (F L ) and the gravitational forces. The theoretical lift forces seem to be only around 10% of the gravitational forces calculated, while they must have been in balance. The relation between the balancing force (F b ) and F G . The analyses revealed a significant force is missing for 5 μm, as big as F G itself. (d) The relation between the balancing force (F b ) and F L .